## How do people evaluate foreign aid to “nasty” regimes?
## Tobias Heinrich & Yoshiharu Kobayashi
## British Journal of Political Science
#########################################################


## Graph for coefficients

## Load estimated models
########################
load("output/ExpModels.Rdata")


## Main graph for A1, different weights
#######################################
out <- vector("list", 3)
for(j in 1:3) out[[j]] <- list(coef(all_exp_mod[[1]][[j]]$Model), sqrt(diag(all_exp_mod[[1]][[j]]$BS_VCV)),
                               rmvnorm(1000, coef(all_exp_mod[[1]][[j]]$Model), all_exp_mod[[1]][[j]]$BS_VCV))
                                             

## Prep the data
coef <- ldply(.data=out, .fun=function(x) data.frame(Name=names(x[[1]]),
                                                     Coefficient=x[[1]],
                                                     Lower=x[[1]] - 1.96 * x[[2]],
                                                     Upper=x[[1]] + 1.96 * x[[2]]))
coef$w <- rep(c(1,2,3), each=nrow(coef)/3)
coef$Name <- rep(c("Intercept", "$25m", "$75m", "Small CT",
                   "Large CT", "Small AML", "Large AML", "Placebo", 
                   "Aid theft", "Rigged election", "Torture", "Media crackdown", 
                   "Aid theft\n   no remedy", "Rigged election\n   no remedy",
                   "Torture\n   noremedy", "Media crackdown\n   noremedy",
                   "Aid theft\n   via U.S. agency",
                   "Aid theft\n   via NGO",
                   "Aid theft\n   via IO",
                   "Rigged election\n   via U.S. agency",
                   "Rigged election\n   via NGO",
                   "Rigged election\n   via IO",
                   "Torture\n   via U.S. agency",
                   "Torture\n   via NGO",
                   "Torture\n   via IO",
                   "Media crackdown\n   via U.S. agency",
                   "Media crackdown\n   via NGO",
                   "Media crackdown\n   via IO"), 3)
                   
coef[grepl(x=coef$Name, pattern="via"), c("Coefficient", "Lower", "Upper")] <- 10*coef[grepl(x=coef$Name, pattern="via"), c("Coefficient", "Lower", "Upper")]

## Insert missing levels
coef <- rbind(coef,
              data.frame(Name=rep(c("$50m", "Baseline benefits", "No issue"), each=3),
                         Coefficient=0, Lower=0, Upper=0,
                         w=rep(c(1,2,3), 3)))

coef$Header <- 0  
coef <- rbind(coef, 
              data.frame(Name=rep(c("Costs", "Benefits", "Potential issues", "Remedy (in $10m)"), each=3),
                         Coefficient=0, Lower=0, Upper=0, 
                         w=rep(c(1,2,3), 4), Header=1))

coef$Name <- as.character(coef$Name)
coef$Name[coef$Header == 0] <- paste0("   ", coef$Name[coef$Header == 0])

coef <- subset(coef, Name != "   Intercept")

## Reorder variable names
##########################
coef$Name2 <- factor(coef$Name, 
                    levels=rev(c("Costs", "   $25m", "   $50m", "   $75m",
                                 "Benefits", "   Baseline benefits", "   Small AML",
                                 "   Large AML", "   Small CT", "   Large CT",
                                 "Potential issues", "   No issue",
                                 "   Placebo", "   Aid theft", "   Rigged election", 
                                 "   Torture", "   Media crackdown",
                                 "   Aid theft\n   no remedy", "   Rigged election\n   no remedy",
                                 "   Torture\n   noremedy", "   Media crackdown\n   noremedy",
                                 "Remedy (in $10m)",
                                 "   Aid theft\n   via U.S. agency", "   Aid theft\n   via NGO", "   Aid theft\n   via IO",
                                 "   Rigged election\n   via U.S. agency", "   Rigged election\n   via NGO", "   Rigged election\n   via IO",         
                                 "   Torture\n   via U.S. agency", "   Torture\n   via NGO", "   Torture\n   via IO",
                                 "   Media crackdown\n   via U.S. agency",
                                 "   Media crackdown\n   via NGO", "   Media crackdown\n   via IO")))
                             
## Names for weighting schemes
coef$weight <- "Complex"
coef$weight[coef$w == 2] <- "Basic + demographics"
coef$weight[coef$w == 3] <- "Basic"


## Graph coefficients for experimental results
g <- ggplot(data=coef, aes(x=Name2, y=Coefficient))
g <- g + geom_point(colour=NA, fill=NA)+ coord_flip() + xlab("") + theme_bw()
g <- g + geom_hline(yintercept=0)
g <- g + facet_grid(~ weight)
g <- g + theme(strip.text.x=element_text(hjust=0, vjust=.2, size=rel(1.23)),
               strip.text.y=element_text(hjust=0, size=rel(1.23)),
               strip.background=element_blank(),
               axis.text = element_text(hjust=0, size=rel(.77)),
               axis.text.y = element_text(hjust=0),
               panel.grid.major.y = element_line(size=0.25, colour="grey80", 
                                                 linetype="dashed"),
               panel.grid.major.x = element_blank(),
               panel.grid.minor = element_blank(),
               axis.ticks = element_blank(),
               plot.title = element_text(size=rel(1.8), hjust=0))
g <- g + geom_linerange(data=subset(coef, Header == 0),
                        aes(x=Name2, ymin=Lower, ymax=Upper),
                        position=position_dodge(width=.8))
g <- g + geom_point(data=subset(coef, Header == 0),
                    aes(x=Name2, y=Coefficient),
                    position=position_dodge(width=.8), size=2.4)
print(g)
ggsave(file="output/figures/ExpRegression-A1-Long.pdf", width=12, height=11)



## Main graph for A2, different weights
#######################################
out <- vector("list", 3)
for(j in 1:3) out[[j]] <- list(coef(all_exp_mod[[2]][[j]]$Model), sqrt(diag(all_exp_mod[[2]][[j]]$BS_VCV)),
                               rmvnorm(1000, coef(all_exp_mod[[2]][[j]]$Model), all_exp_mod[[2]][[j]]$BS_VCV))


## Prep the data
coef <- ldply(.data=out, .fun=function(x) data.frame(Name=names(x[[1]]),
                                                     Coefficient=x[[1]],
                                                     Lower=x[[1]] - 1.96 * x[[2]],
                                                     Upper=x[[1]] + 1.96 * x[[2]]))
coef$w <- rep(c(1,2,3), each=nrow(coef)/3)
coef$Name <- rep(c("Intercept", "$25m", "$75m", "Small CT",
                   "Large CT", "Small AML", "Large AML", "Placebo", 
                   "Aid theft", "Rigged election", "Torture", "Media crackdown"), 3)

## Insert missing levels
coef <- rbind(coef,
              data.frame(Name=rep(c("$50m", "Baseline benefits", "No issue"), each=3),
                         Coefficient=0, Lower=0, Upper=0,
                         w=rep(c(1,2,3), 3)))

coef$Header <- 0  
coef <- rbind(coef, 
              data.frame(Name=rep(c("Costs", "Benefits", "Potential issues"), each=3),
                         Coefficient=0, Lower=0, Upper=0, 
                         w=rep(c(1,2,3), 3), Header=1))

coef$Name <- as.character(coef$Name)
coef$Name[coef$Header == 0] <- paste0("   ", coef$Name[coef$Header == 0])

coef <- subset(coef, Name != "   Intercept")

## Reorder variable names
coef$Name2 <- factor(coef$Name, 
                     levels=rev(c("Costs", "   $25m", "   $50m", "   $75m",
                                  "Benefits", "   Baseline benefits", "   Small AML",
                                  "   Large AML", "   Small CT", "   Large CT",
                                  "Potential issues", "   No issue",
                                  "   Placebo", "   Aid theft", "   Rigged election", 
                                  "   Torture", "   Media crackdown")))

## Names for weighting schemes
coef$weight <- "Complex"
coef$weight[coef$w == 2] <- "Basic + demographics"
coef$weight[coef$w == 3] <- "Basic"


## Graph coefficients for experimental results
g <- ggplot(data=coef, aes(x=Name2, y=Coefficient))
g <- g + geom_point(colour=NA, fill=NA)+ coord_flip() + xlab("") + theme_bw()
g <- g + geom_hline(yintercept=0)
g <- g + facet_grid(~ weight)
g <- g + theme(strip.text.x=element_text(hjust=0, vjust=.2, size=rel(1.23)),
               strip.text.y=element_text(hjust=0, size=rel(1.23)),
               strip.background=element_blank(),
               axis.text = element_text(hjust=0, size=rel(.77)),
               axis.text.y = element_text(hjust=0),
               panel.grid.major.y = element_line(size=0.25, colour="grey80", 
                                                 linetype="dashed"),
               panel.grid.major.x = element_blank(),
               panel.grid.minor = element_blank(),
               axis.ticks = element_blank(),
               plot.title = element_text(size=rel(1.8), hjust=0))
g <- g + geom_linerange(data=subset(coef, Header == 0),
                        aes(x=Name2, ymin=Lower, ymax=Upper),
                        position=position_dodge(width=.8))
g <- g + geom_point(data=subset(coef, Header == 0),
                    aes(x=Name2, y=Coefficient),
                    position=position_dodge(width=.8), size=2.4)
print(g)
ggsave(file="output/figures/ExpRegression-A2-Long.pdf", width=12, height=11)


## Graph that'll go into the manuscript
g <- ggplot(data=subset(coef, weight == "Complex"), aes(x=Name2, y=Coefficient))
g <- g + geom_point(colour=NA, fill=NA)+ coord_flip() + xlab("") + theme_bw()
g <- g + geom_hline(yintercept=0) + ggtitle("Costs, benefits, and potential issues")
g <- g + theme(strip.text.x=element_text(hjust=0, vjust=.2, size=rel(1.23)),
               strip.text.y=element_text(hjust=0, size=rel(1.23)),
               strip.background=element_blank(),
               axis.text = element_text(hjust=0, size=rel(.77)),
               axis.text.y = element_text(hjust=0),
               panel.grid.major.y = element_line(size=0.25, colour="grey80", 
                                                 linetype="dashed"),
               panel.grid.major.x = element_blank(),
               panel.grid.minor = element_blank(),
               axis.ticks = element_blank(),
               plot.title = element_text(size=rel(1.23), hjust=0))
g <- g + geom_linerange(data=subset(coef, Header == 0 & weight == "Complex"),
                        aes(x=Name2, ymin=Lower, ymax=Upper),
                        position=position_dodge(width=.8))
g <- g + geom_point(data=subset(coef, Header == 0 & weight == "Complex"),
                    aes(x=Name2, y=Coefficient),
                    position=position_dodge(width=.8), size=2.4)
ggsave(file="output/figures/ExpRegression-A2-Complex.pdf", width=10, height=5.6)


## Save output
tmp <- subset(coef, weight == "Complex")
tmp[, c("Coefficient", "Lower", "Upper")] <- round(tmp[, c("Coefficient", "Lower", "Upper")], di=2)
write.csv(x=tmp, file=paste("output/tables/A2-Complex Coefficients.csv"))


## Main graph for A3, different weights
#######################################
out <- vector("list", 3)
for(j in 1:3) out[[j]] <- list(coef(all_exp_mod[[3]][[j]]$Model), sqrt(diag(all_exp_mod[[3]][[j]]$BS_VCV)),
                               rmvnorm(1000, coef(all_exp_mod[[3]][[j]]$Model), all_exp_mod[[3]][[j]]$BS_VCV))


## Prep the data
coef <- ldply(.data=out, .fun=function(x) data.frame(Name=names(x[[1]]),
                                                     Coefficient=x[[1]],
                                                     Lower=x[[1]] - 1.96 * x[[2]],
                                                     Upper=x[[1]] + 1.96 * x[[2]]))
coef$w <- rep(c(1,2,3), each=nrow(coef)/3)
coef$Name <- rep(c("Intercept", "$25m", "$75m", "Small CT",
                   "Large CT", "Small AML", "Large AML", "Placebo", 
                   "Aid theft", "Rigged election", "Torture", "Media crackdown",
                   "Aid theft + Small AML", "Rigged election + Small AML", "Torture + Small AML", "Media crackdown + Small AML",
                   "Aid theft + Large AML", "Rigged election + Large AML", "Torture + Large AML", "Media crackdown + Large AML",
                   "Aid theft + Small CT", "Rigged election + Small CT", "Torture + Small CT", "Media crackdown + Small CT",
                   "Aid theft + Large CT", "Rigged election + Large CT", "Torture + Large CT", "Media crackdown + Large CT"), 3)

## Insert missing levels
coef <- rbind(coef,
              data.frame(Name=rep(c("$50m", "Baseline benefits", "No issue"), each=3),
                         Coefficient=0, Lower=0, Upper=0,
                         w=rep(c(1,2,3), 3)))

coef$Header <- 0  
coef <- rbind(coef, 
              data.frame(Name=rep(c("Costs", "Benefits and potential issues"), each=3),
                         Coefficient=0, Lower=0, Upper=0, 
                         w=rep(c(1,2,3), 2), Header=1))

coef$Name <- as.character(coef$Name)
coef$Name[coef$Header == 0] <- paste0("   ", coef$Name[coef$Header == 0])

coef <- subset(coef, Name != "   Intercept")

## Reorder variable names
coef$Name2 <- factor(coef$Name, 
                     levels=rev(c("Costs", "   $25m", "   $50m", "   $75m",
                                  "Benefits and potential issues", "   Baseline benefits", "   Small AML",
                                  "   Large AML", "   Small CT", "   Large CT",
                                  "   No issue",
                                  "   Placebo", "   Aid theft", "   Rigged election", 
                                  "   Torture", "   Media crackdown",
                                  "   Aid theft + Small AML", "   Rigged election + Small AML", "   Torture + Small AML", "   Media crackdown + Small AML",
                                  "   Aid theft + Large AML", "   Rigged election + Large AML", "   Torture + Large AML", "   Media crackdown + Large AML",
                                  "   Aid theft + Small CT", "   Rigged election + Small CT", "   Torture + Small CT", "   Media crackdown + Small CT",
                                  "   Aid theft + Large CT", "   Rigged election + Large CT", "   Torture + Large CT", "   Media crackdown + Large CT")))

## Names for weighting schemes
coef$weight <- "Complex"
coef$weight[coef$w == 2] <- "Basic + demographics"
coef$weight[coef$w == 3] <- "Basic"


## Graph coefficients for experimental results
g <- ggplot(data=coef, aes(x=Name2, y=Coefficient))
g <- g + geom_point(colour=NA, fill=NA)+ coord_flip() + xlab("") + theme_bw()
g <- g + geom_hline(yintercept=0)
g <- g + facet_grid(~ weight)
g <- g + theme(strip.text.x=element_text(hjust=0, vjust=.2, size=rel(1.23)),
               strip.text.y=element_text(hjust=0, size=rel(1.23)),
               strip.background=element_blank(),
               axis.text = element_text(hjust=0, size=rel(.77)),
               axis.text.y = element_text(hjust=0),
               panel.grid.major.y = element_line(size=0.25, colour="grey80", 
                                                 linetype="dashed"),
               panel.grid.major.x = element_blank(),
               panel.grid.minor = element_blank(),
               axis.ticks = element_blank(),
               plot.title = element_text(size=rel(1.8), hjust=0))
g <- g + geom_linerange(data=subset(coef, Header == 0),
                        aes(x=Name2, ymin=Lower, ymax=Upper),
                        position=position_dodge(width=.8))
g <- g + geom_point(data=subset(coef, Header == 0),
                    aes(x=Name2, y=Coefficient),
                    position=position_dodge(width=.8), size=2.4)
ggsave(file="output/figures/ExpRegression-A3-Long.pdf", width=12, height=11)



## Main graph for A4, different weights
#######################################
out <- vector("list", 3)
for(j in 1:3) out[[j]] <- list(coef(all_exp_mod[[4]][[j]]$Model), sqrt(diag(all_exp_mod[[4]][[j]]$BS_VCV)),
                               rmvnorm(1000, coef(all_exp_mod[[4]][[j]]$Model), all_exp_mod[[4]][[j]]$BS_VCV))


## Prep the data
coef <- ldply(.data=out, .fun=function(x) data.frame(Name=names(x[[1]]),
                                                     Coefficient=x[[1]],
                                                     Lower=x[[1]] - 1.96 * x[[2]],
                                                     Upper=x[[1]] + 1.96 * x[[2]]))
coef$w <- rep(c(1,2,3), each=nrow(coef)/3)
coef$Name <- rep(c("Intercept", "$25m", "$75m", "Small CT",
                   "Large CT", "Small AML", "Large AML", "Placebo", 
                   "Aid theft", "Rigged election", "Torture", "Media crackdown",
                   "Aid theft + $25m", "Rigged election + $25m", "Torture + $25m", "Media crackdown + $25m",
                   "Aid theft + $75m", "Rigged election + $75m", "Torture + $75m", "Media crackdown + $75m"), 3)

## Insert missing levels
coef <- rbind(coef,
              data.frame(Name=rep(c("$50m", "Baseline benefits", "No issue",
                                    "Aid theft + $50m", "Rigged election + $50m", "Torture + $50m", "Media crackdown + $50m"), each=3),
                         Coefficient=0, Lower=0, Upper=0,
                         w=rep(c(1,2,3), 7)))

coef$Header <- 0  
coef <- rbind(coef, 
              data.frame(Name=rep(c("Benefits", "Costs and potential issues"), each=3),
                         Coefficient=0, Lower=0, Upper=0, 
                         w=rep(c(1,2,3), 2), Header=1))

coef$Name <- as.character(coef$Name)
coef$Name[coef$Header == 0] <- paste0("   ", coef$Name[coef$Header == 0])

coef <- subset(coef, Name != "   Intercept")

## Reorder variable names
coef$Name2 <- factor(coef$Name, 
                     levels=rev(c("Benefits", "   Baseline benefits", "   Small AML",
                                  "   Large AML", "   Small CT", "   Large CT",
                                  "Costs and potential issues", "   $25m", "   $50m", "   $75m",
                                  "   No issue",
                                  "   Placebo", "   Aid theft", "   Rigged election", 
                                  "   Torture", "   Media crackdown",
                                  "   Aid theft + $25m", "   Aid theft + $50m", "   Aid theft + $75m",
                                  "   Rigged election + $25m", "   Rigged election + $50m", "   Rigged election + $75m", 
                                  "   Media crackdown + $25m", "   Media crackdown + $50m", "   Media crackdown + $75m",
                                  "   Torture + $25m",   "   Torture + $50m",  "   Torture + $75m")))

## Names for weighting schemes
coef$weight <- "Complex"
coef$weight[coef$w == 2] <- "Basic + demographics"
coef$weight[coef$w == 3] <- "Basic"


## Graph coefficients for experimental results
g <- ggplot(data=coef, aes(x=Name2, y=Coefficient))
g <- g + geom_point(colour=NA, fill=NA)+ coord_flip() + xlab("") + theme_bw()
g <- g + geom_hline(yintercept=0)
g <- g + facet_grid(~ weight)
g <- g + theme(strip.text.x=element_text(hjust=0, vjust=.2, size=rel(1.23)),
               strip.text.y=element_text(hjust=0, size=rel(1.23)),
               strip.background=element_blank(),
               axis.text = element_text(hjust=0, size=rel(.77)),
               axis.text.y = element_text(hjust=0),
               panel.grid.major.y = element_line(size=0.25, colour="grey80", 
                                                 linetype="dashed"),
               panel.grid.major.x = element_blank(),
               panel.grid.minor = element_blank(),
               axis.ticks = element_blank(),
               plot.title = element_text(size=rel(1.8), hjust=0))
g <- g + geom_linerange(data=subset(coef, Header == 0),
                        aes(x=Name2, ymin=Lower, ymax=Upper),
                        position=position_dodge(width=.8))
g <- g + geom_point(data=subset(coef, Header == 0),
                    aes(x=Name2, y=Coefficient),
                    position=position_dodge(width=.8), size=2.4)
ggsave(file="output/figures/ExpRegression-A4-Long.pdf", width=12, height=11)



## Calculate differences between gamma_4 and "optimal gamma_5"
##############################################################
out <- vector("list", 3)
for(j in 1:3) out[[j]] <- list(coef(all_exp_mod[[1]][[j]]$Model), sqrt(diag(all_exp_mod[[1]][[j]]$BS_VCV)),
                               rmvnorm(1000, coef(all_exp_mod[[1]][[j]]$Model), all_exp_mod[[1]][[j]]$BS_VCV))
output <- c()

ra_matrix <- matrix(1:25, nrow(out[[1]][[3]]), ncol=25, byrow=TRUE)
for(i in 1:3)
{
  for(j in 1:4)
  {
    for(m in 1:3)
    {
      gamma_4 <- matrix(out[[i]][[3]][, paste0("I(PI == ", j+1, " & RA == 0)TRUE")],
                        nrow=nrow(out[[i]][[3]]), ncol=25)
      gamma_5 <- matrix(out[[i]][[3]][, c(paste0("I(I(PI == ", j+1, " & RT == ", m, ") * RA)"))],
                        nrow=nrow(out[[i]][[3]]), ncol=25)
      diff <- gamma_5 * ra_matrix - gamma_4
      output <- rbind(data.frame(PI=j+1, w=i, Agency=m,
                                 q_m=median(rowMaxs(diff)),
                                 q_l=quantile(rowMaxs(diff), .025), q_u=quantile(rowMaxs(diff), .975)),
                      output)
    }
  }
}

output$w2 <- "Complex"
output$w2[output$w == 2] <- "Basic + demographics"
output$w2[output$w == 3] <- "Basic"
output$PI2 <- "Aid theft"
output$PI2[output$PI == 3] <- "Rigged elections"
output$PI2[output$PI == 4] <- "Torture"
output$PI2[output$PI == 5] <- "Media crackdown"
output$Agency2 <- "U.S. agency"
output$Agency2[output$Agency == 2] <- "NGO"
output$Agency2[output$Agency == 3] <- "IO"

for(i in 1:3)
{
  g <- ggplot(data=subset(output, w == i),
              aes(x=Agency2, ymin=q_l, y=q_m, ymax=q_u))
  g <- g + geom_hline(yintercept=0)
  g <- g + facet_grid(~ PI2) + theme_bw() + coord_flip()
  g <- g + theme(strip.text.x=element_text(hjust=0, vjust=.2, size=rel(1.23)),
                 strip.text.y=element_text(hjust=0, size=rel(1.23)),
                 strip.background=element_blank(),
                 axis.text = element_text(hjust=0, size=rel(.77)),
                 axis.text.y = element_text(hjust=0),
                 panel.grid.major.y = element_line(size=0.25, colour="grey80", 
                                                   linetype="dashed"),
                 panel.grid.major.x = element_blank(),
                 panel.grid.minor = element_blank(),
                 axis.ticks = element_blank(),
                 plot.title = element_text(size=rel(1.8), hjust=0))
  g <- g + geom_linerange(size=.7) + xlab("") + ylab("Difference")
  g <- g + geom_point(size=3)
  ggsave(file=paste0("output/figures/OptAddressbyChannel-w", i, ".pdf"), width=11, height=5)
}  
  
  


